//:Multigrid.h

#ifndef MULTIGRID_H
#define MULTIGRID_H
#include <iostream>
#include <iomanip>
#include <lapacke.h>
#include <cmath>
#include <fstream>
#include <string>

class Restrict{
    void fullweighting(double *f_2h,double *r_h,int m);
    void injection(double *f_2h,double *r_h,int m);
public:
    bool ISFW = true;
    void run(double *f_2h, double *r_h, int m);
};

class Interpolate{
    void linear(double *e_h,double *v_2h,int m);
    void quadratic(double *e_h,double *v_2h,int m);
public:
    bool ISL = true;
    void run(double *e_h,double *v_2h,int m);
};

class Iteration{
public:
    double w = 2.0/3;
    void weightedJacobi(double *v_h,double *f_h,int m);
};

class Bottomsolver{
public:
    void lapack(double *v_h,double *e_h,int m);
};

struct Operator{   
    Restrict res;
    Interpolate inte;
    double maxnorm(double *v_h,int m); 
    void residual(double *r_h,double *f_h,double *v_h,int m);
    void correct(double *v_h,double *e_h,int m);
};

class Multigrid{
    int n,l;
    double **v,**f,*r,*e;   
    Operator op;
    Iteration iteration;
    Bottomsolver bottom;
    void initialize(double *v_h,double *f_h);   
    void VCycle(int m,int num);
    void FullMultigridCycle(int m,int num);
public:
    double V1 = 3,V2 = 3;//times of iteration
    double w = 2.0/3;//weighting factor of weighted Jacobi
    bool ISFW = true,ISL = true,ISVC = true;
    Multigrid(int N); 
    void solve(double *v,double *f,double *u0,int M,double eps,std::ofstream& out);
    ~Multigrid();
};


#endif //MULTIGRID_H
